/* Copyright (C) Nial Peters 2014
 *
 * This file is part of pydoas.
 *
 * pydoas is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * pydoas is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with pydoas.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "fitting.h"


///////////////////////////////////////////////////////////////////////////////

fit_data_t *
createFitDataObj  (double **refs, int num_refs, double **ref_wavelengths,
		           int *ref_lengths, double *fitting_wavelengths,
		           int fit_window_length, int poly_order){

	register int i,j; //counters
	double *first_order, *current_order, *previous_order;
	fit_data_t *obj = NULL;

	if ((poly_order < 0) || (poly_order > 20)){
		//20 is arbitrary - but you shouldn't be fitting with polynomials
		//of this high order.
		ERROR(ERR_VAL, "Background polynomial must have an order in the range 0-20.");
		return NULL;
	}

	if ((obj = malloc(sizeof(fit_data_t))) == NULL){
		ERROR(ERR_MEM, "Failed to allocate memory for fit data object.");
		return NULL;
	}

	//copy the values we need to keep track of
	obj->poly_order = poly_order;

	//create the shift and squeeze object
	obj->sns_data = createSNSDataObj(refs, num_refs, ref_wavelengths,
			                         ref_lengths, fitting_wavelengths,
			                         fit_window_length);

	if (obj->sns_data == NULL){
		ERROR(ERR_MISC, "Failed to create the shift and squeeze data object.");
		destroyFitDataObj(&obj);
		return NULL;
	}

	//allocate the memory for the fitting matrix. Need to +1 to poly_order
	//to account for 0'th order part (constant offset)
	obj->fit_matrix = malloc((num_refs + poly_order + 1) *
			                 fit_window_length * sizeof(double));

	if (obj->fit_matrix == NULL){
		ERROR(ERR_MEM, "Failed to allocate memory for fitting matrix.");
		destroyFitDataObj(&obj);
		return NULL;
	}

	//now populate the polynomial part of the fit_matrix - the rest is
	//filled in when a fit is done
	i = 0;
	//0th order
	for (i=0; i < fit_window_length; i++){
		obj->fit_matrix[i] = 1.0;
	}

	//1st order
	first_order = obj->fit_matrix+fit_window_length;
	if (poly_order > 0){
		memcpy(first_order, fitting_wavelengths, fit_window_length * sizeof(double));
	}

	//subsequent polynomial orders
	for (j=2; j<=poly_order; j++){
		current_order = obj->fit_matrix + (j*fit_window_length);
		previous_order = obj->fit_matrix + ((j-1)*fit_window_length));

		for (i=0; i<fit_window_length; i++){

			current_order[i] = previous_order[i] * first_order[i];
		}
	}

	return obj;

}

///////////////////////////////////////////////////////////////////////////////

void
destroyFitDataObj (fit_data_t **data){
	sns_data_t *d = *data;

	if (d == NULL){
		return;
	}

	destroySNSDataObj(&(d->sns_data));

	if (d->fit_matrix != NULL){
		free(d->fit_matrix);
		d->fit_matrix = NULL;
	}

	//finally set the pointer to NULL so that it cannot be used by mistake
	*data = NULL;
}

///////////////////////////////////////////////////////////////////////////////

static inline void
fittingFunction(double *p, double *x, int m, int n, void *data){
	/*
	 *
	 */
	register int i;
	fit_data_t *fd = (fit_data_t *)data;
	const int num_refs = fd->sns_data->num_refs;
	const int fit_window_length = fd->sns_data->fit_window_length;
	double chisq;

	//shift and squeeze the references by the desired amounts
	for(i=0; i<num_refs; i++){
		shift_and_squeeze(fd->sns_data, p[i*2], p[(i*2) + 1], i, fd->ref_ptrs[i]);
	}

    //TODO - everything below here needs to be finished
	//gsl_vector_view fit_result_vec = gsl_vector_view_array(x, fit_window_length);

	//gsl_multifit_linear(&(fd->lin_fit_params).matrix, &fit_result_vec.vector,fd->lin_fit_result, fd->cov, &chisq, fd->workspace);

	//gsl_blas_dgemv(CblasNoTrans, 1.0, &(fd->lin_fit_params).matrix, fd->lin_fit_result, -1.0, &fit_result_vec.vector);

	//result is now stored in x, as required by levmar fitting function
}

///////////////////////////////////////////////////////////////////////////////

int
fitSpectrum(fit_data_t *fit_obj, double *wavelengths, double *counts,
		    int spec_length){




	return 0; //TODO
}
